Density functional theory description of hole-trapping in Si02: a successful 

self-interaction-corrected approach. 
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We present a self-interaction-corrected (SIC) density-functional-theory (DFT) approach for the 
description of systems with an unpaired electron or hole such as spin 1/2 defect-centers in solids 
or radicals. Our functional is easy-to-implement and its minimization does not require additional 
computational effort with respect to ordinary DFT functionals. In particular it does not present 
multi-minima, as the conventional SIC functionals. We successfully validate the method studying 
the hole self-trapping in quartz associated to the Al substitutional impurity. We show that our 
approach corrects for the well known failures of standard DFT functionals in this system. 

PACS numbers: 71.15.Mb, 61.72.Bb, 71.55-i 
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When an electron, a hole, or an electron-hole pair is 
introduced in an insulator, its wavefunction can have dif- 
ferent degrees of localization. It can be localized (self- 
trapped) on a single atom/bond or it can be delocalized 
on a larger scale. Trapped electron or hole centers oc- 
cur in very different systems. The F-centers in alkali 
halides are paramagnetic centers formed by electron de- 
fects trapped on negative ion vacancies in the crystal. 
They can be obtained by introducing a stoichiometric 
excess of alkali metal atoms or ionizing the crystal with 
radiationP|. Trapped hole defects are found in alkaline 
earths oxides ([M]° centers) "i^ upon substitution of the 
divalent alkaline earth atom with a monovalent one. 

Several examples of self-trapping occur in silica. Holes 
and electrons induced by ionizing radiation self-trap in 
amorphous 8102 01 but not in quartzj^. In quartz, lo- 
calization of holes is achieved via the substitution of Si 
with a trivalent atom (e.g. Al)|^. Similarly, electron 
self-trapping is obtained substituting a Si with a pen- 
tavalent atom (e.g. P). Self-trapped holes or electrons in 
quartz can also be obtained by substituting Si with Ge 
and by removing or adding an electron with ionizing ra- 
diation. Excitons self-trap both in amorphous Silica and 
in quartz 

The study of self-trapping of centers in both amor- 
phous and crystalline Si02 is a subject of important tech- 
nological implications. The creation of self-trapped cen- 
ters in Si02 due to ionizing radiation or high intensity 
UV light determines the degradation of UV-transmitting 
fibers(for e.g. UV lithography). Moreover self-trapped 
defects are partly responsibles for the failure rate of MOS 
devices which use amorphous Si02 as an insulating layer. 

From a theoretical point of view, the description of 
trapped defects is particularly challenging since a stan- 
dard density functional theory (DFT) approach based on 
the local spin density approximations (LSDA) or on its 
improved version, the spin polarized generalized gradi- 
ents approximations (SPGGA), often fails to reproduce 
the localization of the defect wavefunction. Unrestricted 
Hartree Fock (UHF) calculations 0, 01 usually correctly 



reproduces the self-trapping. However the UHF descrip- 
tion of the defect-free system is usually less accurate than 
that obtained with DFT. Also, UHF is computationally 
more demanding. 

DFT predicts most F-centers in alkali halides to be 
delocalized ..ij. In the case of Si02, electron self-trapping 
is correctly described by DFT, whereas holes are found 
delocalized 0, A prototype example are neutral Al 
centers in Si02. The structure of pure silica is composed 
of corner-sharing [Si04] tetrahedrons. The Al doping oc- 
curs as a substitution of a Si in the tetrahedral structure, 
resulting in consequential modifications of electronic and 
geometric properties. The Al-Si substitution introduces 
in the system a hole which, according to electron spin 
resonance, is localized on one of the surrounding oxygen. 
On the contrary, and in clear disagreement with experi- 
mental findings, DFT using standard functionals predicts 
a hole wavefunction delocalized over the four surround- 
ing oxygens 6, 8J. For this reason the description of hole 
trapping on Al impurities in Silica has been defined "a 
challenge for density functional theories"Ja| • 

As suggested by several authors 0, Q the failure of 
DFT in describing self-trapping in general, and self- 
trapping on Al centers in particular, might be due to 
the incomplete cancellation of the unpaired-electron self- 
interaction. This cancellation occurs exactly in UHF. A 
possible solution to the problem might be to use self- 
interaction corrected (SIC) functionals^, It is in- 
deed well known that SIC functionals can describe elec- 
tronic states which are not reproduced by LSDA/SPGGA 
functionals jllj. Unfortunately, the implementation of 
DFT functionals with self-interaction corrections on each 
orbital, besides being technically complex, leads to multi- 
minima problems and sometimes degrades the results ob- 
tained with standard functionals |l2l |. 

In this work we show that SIC functionals are indeed a 
remedy to the failure of DFT in describing Al defects in 
Silica. We present an easy-to-implement SIC functional 
for one-particle spin 1/2 defects in solids. In our SIC 
functional approach we apply the self-interaction correc- 
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tion to the unpaired electron only. We treat the closed 
shell system formed by the remaining 2N electrons using 
standard DFT functionals but imposing the conditions 
that up and down electrons with the same orbital quan- 
tum numbers have identical spatial wavefunctions. Our 
method is free from the multi-minima problem found in 
previous SIC implementations and its computational cost 
is equal to that of ordinary DFT functionals. We validate 
the method against experimental data by calculating op- 
timized geometries and hyperfine coupling parameters of 
Al defects in silica. 

We consider a system oi 2N + 1 electrons with Nt^ — 
-I- 1 up electrons and = N down electrons. The 
wavefunction of the i*'' a— electron (a =t, J.) is \tpicr)- 
Throughout the paper we assume the wavefunctions \ipia) 
to be orthonormal. We write ricr (r) = X^ilTi \ and 
'^(i") = Sdr'^cr(r). The magnetization density is m(r) = 
n^ir) -ni{r). 

In spin polarized density functional theory the total 
energy functional is: 

F[{^P,„}] = T[{^j,,}]+EKs[n^,n^] + 

+ E E Kj^al^ja) - S^,) (1) 

where T[{^^a}] = -^EaEf' is the sin- 
gle particle kinetic energy and EKs[ni,ni] = Ecxt[n] + 
E-aln] -I- £^xc['T'T) ^i]- The functionals Ecxt[n], E}i[n] and 
Exc[n-\ ,ni] are the external potential and the Hartree 
and exchange-correlation functionals respectively 0. 
Atomic units are used throughout the paper. Lagrange 
multipliers A°!^ are included to impose the wavefunction 
orthonormalization condition. The ground state energy 
is obtained by minimizing the total energy functional, 
namely by imposing '^''^^ff '°|"^^ = 0, for each i and a. 
Upon derivation of the total energy functional respect 
to (f/'io-l, we obtain: 
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(5-Eks["t."-i] 



where V^sir) - ^„^(,) ■ 

The Hartree functional En in ^[{7/1^0.}] contains a self- 
interaction term for each electron. If the exact functional 
is used, these terms are exactly canceled by an identical 
term (with opposite sign) included in E^c- When approx- 
imate forms of E^c are used, a self-interaction is intro- 
duced in F[{ipic,}], since E^c no longer cancels the self 
interaction terms in En exactly. This self-interaction is 
unphysical and must be subtracted. Since there is one 
spurious self-interaction term per electron, the elimina- 
tion of all the terms leads to an orbital-dependent correc- 
tion in the functional. The resulting Kohn-Sham equa- 
tions are not invariant anymore for a unitary transforma- 



tion in the subspace of occupied orbitals, and the min- 
imization of the functional leads to multi-minima prob- 
lems. When DFT calculations give a good description 
of the system without the spin 1/2 center (i.e. pure 
8102) the relevant self- interaction is typically the un- 
paired electron's one0, 0. As a consequence we pro- 
pose to subtract the self-interaction only for the un- 
paired electron. More specifically, we remove only the 
self-interaction term associated with the magnetization 
density m(r). Indeed, using standard DFT functionals, 
the Kohn-Sham eigenstates are such that {ipii) — \ipii) 
and m(r) ~ KrjV'Ar+ii)!^ (for which eigenstates we have 
A°j = Sijeia-, where eia- are the Kohn and Sham eigenval- 
ues). As it will be shown below, this strategy eliminates 
both the orbital dependence of the equations and the 
multi-minima problem. 

A straightforward way to subtract, at least partially, 
the self-interaction was developed by Perdew and Zunger 
(SPZ) 1^. Applying this approach to the unpaired- 
electron leads to the definition of a new self-interaction 
corrected functional (denoted SPZ), Fsic'iii^icr}] = 
^[{V'o-}] + ^Fspz[n^,ni], where 



AFspz[nT,ni] = -E'hM - Exc[m,0] 



(3) 



This amounts to subtracting the self-interaction term as- 
sociated with the magnetization density from both the 
Hartree and from the exchange correlation functionals. 

Besides the SPZ scheme, in this work we propose a new 
SIC functional (denoted US), defined as Fsic[{tpia}] = 
F[{il;,^}] + AFus, where 

AF(7s[n|,n|] = -En[m] - E^dn^^ni] + E^dn-^ ~m,ni] 

(4) 

namely, (i) we subtract the unpaired electron self- 
interaction from the Hartree functional and (ii) we re- 
place the exchange correlation functional for the 2N-I-1 
electrons system with the one for the 2N electrons system 
without the unpaired electron. 

Unfortunately, the US or PZ self-interaction corrected 
functionals are not sufficient by themselves to obtain 
physically relevant densities. Due to the quadratic de- 
pendence (with a negative sign) of En respect to m(r), 
the two functionals tend to maximize everywhere |m-(r)| 
by separating the spin up and spin down densities. As 
a consequence becomes very different from 

for i < N and m(r) is not approximately equal to 
|(r|?/'Ar+i|)p. This unphysical solution can be ehminated 
by introducing a second constraint on the 2N-electrons 
system, i.e. the system without unpaired electron. We 
impose that up and down electrons with the same or- 
bital quantum numbers have identical spatial wavefunc- 
tion (spin-restricted solution): 

I^U) for Z=1,...,A^ (5) 

IV'w+it) = \fpN+i) (6) 
In this way the total energy becomes a function of {ipi}: 
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where the 'qij Lagrange multipHers are used to enforce 
the orthonormalization conditions. Ai^sic["-Ti "-il 
SIC correction in the SPZ (eq. ^ or in the US (eq. \^ 
scheme. Minimization of the total energy functional can 
be achieved by imposing — for i — 1,...,N + 1. 

Functional derivation respect to leads to: 



Functional Delocalized geometry Localized geometry 

PBE stable unstable (43.9 mRy) 

SPZ stable metastable (1.3 mRy) 

US unstable (238 mRy) stable 

UHF unknown stable 



TABLE I: Stability of of localized and delocalized geome- 
tries with different functionals. In parenthesis we indicate the 
energy difference between the metastable/unstable structure 
and the stable one. With PBE we use the US minimal energy 
structure as the localized geometry. With US we use the PBE 
minimal energy structure as the delocalized geometry. 
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l^w+i) + (8) 

N+l 

(9) 



The potential Vg^^ is is defined as VgJc(r) — V-^g{r) + 
^st^ll) , where AFsic is Ai^spz or AFus depending on 
which subtraction scheme is adopted (eq. |31or eq. 

Equations |H1 and El form a set of self-consistent equa- 
tions coupled via the terms involving Lagrange multipli- 
ers. The condition of = is equivalent to the mini- 
mization of the corresponding functional. For fixed ionic 
positions we minimize the SPZ and US functionals with 
the constraints in eq. I5l6l using the the grad ients given in 
eq. IHlandlHland the Car-Parrinello|14l Il5 1 m ethod with 



a damped molecular dynamics approach ^ . We used 
combined electronic and ionic dynamics for the geome- 
try optimizations. 

We simulate the structure of an Al defect in Silica us- 
ing a neutral cell of 72 atoms (one Al atom, 23 Si and 48 
O). We perform electronic structure calculations [isllr^ 
using DFT in the spin polarized generalized gradient 
approximation and the Perdew-Burke-Ernzerhof (PBE) 
functional 0| corrected for self-interaction as in eq. El 
and in eq. ^ We use norm conserving pseudo-potentials 
|l9l |. The wave functions are expanded in plane waves 
using a 70 Ry cutoff. We sample the Brillouin zone with 
the P point and we impose the (rjV'i) to be real [20j . 

The results of geometrical optimization using different 
functionals are illustrated in tables HI andllT] We label de- 
localized geometry a structure in which the hole wavefunc- 
tion is delocalized on the four O atoms surrounding the 
Al one. On the contrary a localized geometry is a struc- 
ture in which the hole wavefunction is localized on one of 
the four surrounding Oxygens. For a given functional we 
indicate the structure corresponding to an absolute min- 
imum with stable and to a local minimum (but not the 
absolute one) with metastable. The structure not corre- 
sponding to any minimum is called unstable. The PBE 





Delocalized 


geometries 


Localized geometries 


Bond 


PBEfia] 


SPZ 


SPZ 


US 


UHFf6] 


Al-O(l) 


1.744 


1.742 


1.938 


1.957 


1.924 


Al-0(2) 


1.744 


1.742 


1.713 


1.710 


1.688 


Al-0(3) 


1.753 


1.748 


1.706 


1.705 


1.703 


Al-0(4) 


1.753 


1.748 


1.719 


1.715 


1.689 


Si-O(l) 


1.613 


1.596 


1.756 


1.793 




Si-0(2) 


1.613 


1.596 


1.604 


1.602 




Si-0(3) 


1.618 


1.595 


1.595 


1.595 




Si-0(4) 


1.618 


1.595 


1.601 


1.601 





TABLE II: Bond-lengths (A) around the substitutional Al im- 
purity using unrestricted Hartree-Fock and density functional 
theory with functionals PBEIi3,SPZ and US. 



stable structure j21 is a delocalized geometry, in qualita- 
tive disagreement with UHF. The SPZ stable structure 
is delocalized. However a metastable localized solution 
occurs at slightly higher energy (1.3 mRyd). The US 
stable structure is localized, in agreement with the UHF 
results. We did not find any metastable solution using 
the US functional and in particular the delocalized struc- 
tures found with PBE or SPZ are unstable. Finally, to 
judge the effects of the spin-restricted prescription, we 
compute the total energy difference between the local- 
ized and the delocalized structures using PBE with the 
constraint of eqs. [31 and We obtain 45.6 mRyd which 
is very close to the unrestricted result of 43.9 mRyd given 
in tabled i.e. the spin-restricted condition weakly affects 
the total energy differences. 

Delocalized geometries are tetrahedral structures with 
Al-0 and Si-0 bond-lengths very close to the Al-0 dis- 
tance (« I.73A) in AIPO4 and to the Si-0 bond-lengths 
in quartz (« 1.61 A), respectively. The localization of the 
hole on one particular oxygen (labeled 0(1) in tab. ^ 
leads to a distorted tetrahedral structure with two elon- 
gated Al-O(l) and Si-O(l) bonds, while all the others 
Al-0 bonds are only slightly smaller than in the delo- 
calized case. The Al-0 and Si-0 bond-lengths of stable 
US and metastable SPZ structures are in good agreement 
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with UHF. 



Geom; 


PBE 


US 


SPZ 


US 


UHFf6] 




Fun: 


PBE 


PBE 


SPZ 


US 


UHF[6] 


Exp. 




-17.0 


-72.9 


-96.9 


-89.5 


-89.3 


-85.0 


0(1) 


8.46 


36.3 


48.3 


44.6 


44.6 


41.2 




8.55 


36.7 


48.6 


44.8 


44.7 


43.8 



TABLE III: Anisotropic hyperfine parameters (Gauss) of the 
^^O (1) atom. The first row "Geom:" indicates the functional 
used to determine the geometry, the second row "Fun:" the 
functional used in the hyperfine coupling calculation. In this 
table SPZ geometry refers to the metastable localized struc- 
ture obtained with the SPZ functional (see text). 

To validate our SIC methods against experiment we 
compute hyperfine parameters using the formalism devel- 
oped in ref. |2]|. The hyperfine interaction includes an 
isotropic part (Fermi contact) and an anisotropic (dipo- 
lar) part. The isotropic part measures the spin density 
at the nucleus and it is non-zero only for wavefunctions 
containing s-wave components [U. The hole state, both 
in the localized and delocalized solutions, corresponds 
to the O 2p lone-pair state. Thus the Fermi contact is 
mainly due to the spin polarization of the doubly oc- 
cupied O 2s state. Since the spin-restricted conditions 
(eqs. |5l and 13 suppresses such spin polarization, the 
Fermi contact term cannot be computed with our spin- 
restricted SIC functionals On the other hand, we 
can access the anisotropic terms since they capture the 
p— like component of the electron wavefunction which is 
weakly affected by the spin-restricted condition In 
table lllll we report the principal values of the dipolar part 
for ^'^O atoms. Our implementation of the SIC functional 
substantially improves the PBE results, giving a dipolar 
part very close to the UHF results and the experimental 
data. 

In this work we have presented a very effective and 
easy-to-implement self-interaction corrected approach. It 
can be applied to hole or electron spin-l/2-centers in 
solids. We have validate our method by calculating en- 
ergetics, geometries and hyperfine couplings of neutral 
Al substitutional defects in quartz. Our SIC approach 
corrects the known deficiencies 0, |^ of standard DFT 
functionals and gives results in good agreement with ex- 
periments and with the self-interaction free UHF calcu- 
lations. Thus the functional proposed in this work solves 
the problem^ of describing the behavior of hole self- 
trapping in silica in the framework of density functional 
theory. Finally we note that our approach can also be 
applied to molecular systems such as spin 1/2 radicals 
which represents a fundamental subject of research in 
chemistry and biochemistry. 
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